Deterministic Epidemic SEIR Model Fit Using ABC

Introduction

SIR/SEIR models are common epidemic modeling tools, with demonstrated success. The model is a non-linear differential equation given by.

\[ \frac{dS}{dt} = - \beta * S * I / N\] \[ \frac{dE}{dt} = \beta * S * I / N - \sigma * E\] \[ \frac{dI}{dt} = \sigma * E - \gamma * I \] \[ \frac{dR}{dt} = \gamma * I\]

This model uniquely defines an epidemic curve, in fact we can simualte from the curve using the deSolve package in R to get a sense of what it looks like.

library(deSolve)

## Create an SIR function
seir <- function(time, state, parameters) {
  
  with(as.list(c(state, parameters)), {
    
    dS <- -beta * S * I
    dE <-  beta * S * I - sigma * E
    dI <- sigma*E - gamma*I
    dR <-   gamma * I
    
    return(list(c(dS,dE, dI, dR)))
  })
}

### Set parameters
## Proportion in each compartment: Susceptible 0.999999, Infected 0.000001, Recovered 0
init       <- c(S = 1-1e-6, E=0,I = 1e-6, R = 0.0)
## beta: infection parameter; gamma: recovery parameter
parameters <- c(beta = .95, gamma = 1/3.5,sigma=1/3.7)
## Time frame
times      <- seq(0, 100, by = 1)

## Solve using ode (General Solver for Ordinary Differential Equations)
out <- ode(y = init, times = times, func = seir, parms = parameters)
## change to data frame
out <- as.data.frame(out)
plot(out$I,type='l',ylab='Percent Infected',xlab='Time in Days')

We can also think of this model as describing flow between compartments, where we parameterize the rate of flow. The rate of flow between the susceptible and exposed compartment is parameterized by \(\beta\), and the rate of flow between the exposed compartment and the infected compartment by \(\sigma\), and finally the rate of flow between the infected and recovered compartment by \(\gamma\). This is graphically demonstrated below.

SEIR Comparmental Model

This model defines a smooth epidemic curve, which we rarely observe in practice. For instance, here is the current state of the COVID-19 epidemic in NY from https://covidtracking.com/.

## Please visit openintro.org for free statistics materials

## 
## Attaching package: 'openintro'

## The following object is masked from 'package:ggplot2':
## 
##     diamonds

## The following objects are masked from 'package:datasets':
## 
##     cars, trees

location_idx <- 1
for (location in head(locations,51)){
  p <- ggplots[[location_idx]] + labs(title = census_2010[census_2010$ab == location,]$state)
  suppressWarnings(print(p))
  location_idx <- location_idx + 1
}